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Abstract 


This report describes the analysis component of the Goddard Earth Observing Sys- 
tem, Data Assimilation System, Version 1 (GEOS— 1 DAS). The general features of the 
data assimilation system are outlined, followed by a thorough description of the sta- 
tistical interpolation algorithm, including specification of error covariances and quality 
control of observations. We conclude with a discussion of the current status of develop- 
ment of the GEOS data assimilation system. 

The main components of GEOS-1 DAS are an atmospheric general circulation model 
and an Optimal Interpolation algorithm. The system is cycled using the Incremental 
Analysis Update (IAU) technique in which analysis increments are introduced as time 
independent forcing terms in a forecast model integration. The system is capable of pro- 
ducing dynamically balanced states without the explicit use of initialization, as well as a 
time-continuous representation of non-observables such as precipitation and radiational 
fluxes. This version of the data assimilation system was used in the five-year reanalysis 
project completed in April 1994 by Goddard’s Data Assimilation Office (Schubert et 
al. 1993). Data from this reanalysis are available from the Goddard Distributed Ac- 
tive Archive Center (DAAC) which is part of NASA’s Earth Observing System Data 
and Information System (EOSDIS). For information on how to obtain these data sets 
contact the Goddard DAAC at (301) 286-3209, EMAIL daac@gsfc.nasa.gov, or consult 
DAO s Mosaic Home Page URL http://hera.gsfc.nasa.gov/dao.hoine-pag6.htnil. 

This document is available electronically via anonymous ftp from dao . gsf c . nasa . gov, 
directory pub/tech_memos, file volume_4.ps . Z. 
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1 Introduction 


A recent NRC panel report (National Research Council 1991) provides an overview of the 
issues involved in using data assimilation in climate research. In its report the panel recom- 
mends the routine generation of research-quality, model-assimilated and tested geophysical 
data sets to serve a broad range of national endeavors, including climate and global change 
research and prediction. In early 1991 the Data Assimilation Office (DAO) was formed 
within the NASA/Goddard Laboratory for Atmospheres (GLA) and charged with produc- 
ing research-quality analyses for general Earth Science applications. 

In April 1994 the DAO completed a five-year reanalysis for the years 1985 through 1989 
(Schubert et al. 1993). While providing an important tool for climate research, this analysis 
will also serve as a baseline for further assimilation system development. The DAO’s longer 
term goal is to develop an assimilation system capable of making optimal use of the large 
volume of high-quality observations expected near the turn of the century from NASA’s 
Earth Observing System (EOS) satellites. 

The main objective of this Technical Memorandum is to describe the analysis component 
of the system used for the five-year reanalysis mentioned above. This system has been 
named the Goddard EOS Version 1, Data Assimilation System, or GEOS-1 DAS. Its im- 
mediate predecessor is the optimum interpolation (01) analysis scheme described in Baker 
et al. (1987), hereafter denoted B87. The atmospheric general circulation model (AGCM) 
component of the system is referred to as GEOS-1 GCM, and is documented in Volume 1 
of this Technical Memorandum series (Takacs et al. 1994). 

This report describes the analysis component of GEOS-1 DAS with particular emphasis on 
those features which differ significantly from B87. Sections 2-3 present a general outline 
of the data assimilation system and forecast model, followed by a description of the sta- 
tistical interpolation algorithm (section 4), quality control of observations (section 5) and 
localization of the analysis problem (section 6). The specification of forecast and observa- 
tional error statistics is covered in sections 7 and 8. We conclude with a brief discussion 
of the improvements made to GEOS DAS since its version 1. The appendix gives detailed 
derivations of the statistics used in the GEOS-1 DAS and provides a list of acronyms. 

This document is available electronically via anonymous ftp from dao.gsfc.nasa.gov, 
directory pub/tech_memos, file volume_4 .ps .Z or from DAO’s Mosaic Home Page URL 
http : //hera . gsf c . nasa . gov/dao . home .page . html. 


2 General Features of GEOS-1 DAS 

The GEOS-1 DAS analyzes global sea level pressure and near surface winds over the oceans, 
as well as geopotential height, vector wind, and water vapor mixing ratio on constant 
pressure surfaces. The upper air height/wind analyses and the sea level pressure/surface 
wind analyses are done using multivariate statistical interpolation algorithms in which mass 
(height/pressure) and wind data affect both the mass and wind analyses. The moisture 
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analysis is done with a univariate statistical algorithm, and only at levels from 1000 hPa to 
300 hPa. The basic GEOS-1 DAS configuration consists of a 2° latitude by 2.5° longitude, 
14-level analysis (20, 30, 50, 70, 100, 150, 200, 250, 300, 400, 500, 700, 850, 1000 hPa) 
coupled to a 20-level, 2° by 2.5° model for the troposphere and lower stratosphere. Another 
version of the system (not used in the 5-year reanalysis) consists of a 2° by 2.5°, 18-level 
(surface to 0.4 hPa) analysis coupled to a 46-level, T by 2.5° model. This later system has 
been named STRATAN for stratospheric analysis (Coy et al. 1994 and references therein). 
Figure 1 shows the vertical extent and distribution of layers for both the analysis and model 
components of GEOS-1 DAS. 



Figure 1: Vertical extent and distribution of layers for both the analysis and model com- 
ponents of GEOS-1 DAS. Solid lines are analysis levels and hatched lines are model sigma 
levels. 


The analysis scheme 

The GEOS-1 DAS statistical analysis scheme uses up to 75 observations to analyze all 
grid-points within a small three-dimensional cluster. The data selection algorithm, which 
chooses the observations to be used from those that have passed the quality control proce- 
dures, is an empirically tuned decision tree which uses a priori observation error estimates 
in making its choices. All observations used in the analysis have passed a two stage quality 
check (Seablom et al. 1991). The first stage gross check makes use of the assimilation’s 
forecast error variance fields to estimate expected innovation vector variances. In a second 
stage a buddy check compares suspect observations with neighboring data by means of a 
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successive correction analysis to the location of the suspect datum. Forecast error correla- 
tions are modeled with the damped cosine function as in B87, but the fit parameters have 
been recalculated using GEOS-1 DAS. Details are given in section 8 and in the Appendix. 


Schematic of Incremental Analysis Update (IAU) 



Figure 2: Schematic diagram of Incremental Analysis Update (IAU) procedure. 


Incremental Analysis Updates 


The assimilation scheme does not include an explicit initialization step but rather relies 
on the damping properties of a Matsuno time differencing scheme to control initial imbal- 
ances generated by the insertion of analysis increments. However, the initial imbalances 
and spin-up have been greatly reduced over earlier versions by the introduction of an incre- 
mental analysis update (IAU) procedure (Bloom et al. 1991). As shown in Fig. 2, in the 
IAU procedure the analysis increments are computed in a conventional way at the analysis 
times (00, 06, 12, 18 UTC). These increments are then inserted gradually into the model 
by restarting the short-term forecast that provided the analysis background and adding a 
fraction of the analysis increment at each model step. Over the six hour period centered on 
the analysis time, the full effect of the analysis increment is realized. The final assimilation 
product thus effectively consists of a model forecast (heavy solid line in Fig. 2) produced 
using additional heat, momentum, moisture and mass tendency terms which are updated 
every six hours from observations. This update scheme for the assimilation is similar to 
the way the forcing tendencies from the model’s physical parameterizations and filters are 
recomputed intermittently (usually every 30 minutes) and gradually introduced into the 
ongoing integration at every time-step. An important difference between the IAU scheme 
and the usual Newtonian nudging procedure (Anthes 1974; Stauffer and Seaman 1990) is 
that the IAU forcing terms are held constant over the insertion period, while in Newto- 
nian nudging they are proportional to the changing difference between a target analysis 
and the instantaneous current model state. Significant improvements in terms of assimila- 
tion accuracy, noise control, and the hydrologic cycle spin-up are obtained using the IAU 
technique. 
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The actual implementation of IAU in GEOS-1 DAS involves the following steps: 

1. The analysis fields must be interpolated from mandatory levels to the sigma-levels 
required by the GEOS GCM. To minimize the adverse effects of the vertical inter- 
polation, only the analysis increments are vertically interpolated, and then added to 
the original model first guess in sigma-coordinates. For details of the interpolation 
algorithm refer to Takacs et al. (1994). 

2. The analysis variables must be converted to the GCM prognostic variables. For exam- 
ple, the potential temperature field is first computed from analyzed heights and mixing 
ratio. Because the model thermodynamic equation is formulated in flux form (Takacs 
et al. 1994, eq. 9), the after analysis potential temperature must be mass- weighted, 
i.e., the potential temperature field must be multiplied by the factor II = p 3 — pr, 
where p s is the surface pressure and pr = 10 hPa is the pressure at the top of the 
model. 

3. The IAU forcing terms are formed by subtracting the after analysis prognostic fields 
from the corresponding first guess fields, and then dividing this difference by the 
number of seconds in 6 hours. IAU forcing fields are produced for surface pressure, 
wind components, mass-weighted temperature and mass-weighted moisture. 

4. The model integration is restarted 3 hours before the synoptic time and continues for 
6 hours with the IAU forcing held constant in time. After these 6 hours, the forcing 
terms are set to zero and the model integration continues for another 3 hours up to 
next synoptic time when the model state is used as first guess for the next analysis. 


3 The GEOS-1 Model 

This section presents a summary of the main features of the GEOS-1 GCM. This model is 
fully documented by Takacs et al. (1994) which is available electronically via anonymous ftp 
from dao.gsfc.nasa.gov, directory pub/tech-memos, file volume_l.ps. Z or from DAO’s 
Mosaic Home Page URL http://hera.gsfc. nasa.gov/dao. home-page. html. 

The GEOS-1 GCM uses the second-order potential enstrophy and energy conserving hor- 
izontal differencing scheme on a C-grid developed by Sadourney (1975), and further de- 
scribed by Burridge and Haseler (1977). An 8th-order Shapiro filter with a reduced coeffi- 
cient is applied to the wind, potential temperature and specific humidity to avoid non-linear 
computational instability. The reduced-coefficient filter is applied at every step in such a 
way that the amplitude of the two-grid interval wave is essentially removed in six hours. 
Applying the filter weakly at each step way eliminates the shock that occurred in earlier 
assimilations using an intermittent application of the filter. The model also uses a polar 
Fourier filter to avoid linear instability due to violation of the CFL condition for the Lamb 
wave and interval gravity waves. This polar filter is applied only to the tendencies of the 
winds, potential temperature, specific humidity and surface pressure. The model’s vertical 
finite differencing scheme is that of Arakawa and Suarez (1983). The dynamics routines are 
organized into a plug-compatible module called the “ARIES/GEOS dynamical core” which 
is described in Suarez and Takacs (1995). 
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The infrared and solar radiation parameterizations follow closely those described by Harsh- 
vardhan et al. (1987). In the longwave, water vapor absorption is parameterized as in Chou 
(1984), the 15 micron band of C0 2 as in Chou and Peng (1983), and ozone absorption as 
in Rogers (1968) with modifications suggested by Rosenfield et a 1. (1987). The shortwave 
follows Davies (1982), as described in Harsh vardhan et al. (1987). Shortwave absorption by 
water vapor uses a distribution approach as in Lacis and Hansen (1974). Cloud albedo 
and transmissivity for the model layers are obtained from specified single-scattering albedo 
and cloud optical thickness using the delta-Eddington approximation (Joseph et al. 1976; 
King and Harshvardhan 1986). 

The penetrative convection originating in the boundary layer is parameterized using the 
Relaxed Arakawa-Schubert (RAS) scheme (Moorthi and Suarez 1992) which is a simpli- 
fied Arakawa-Schubert (1974) type scheme. As an approximation to the full interaction 
between the different allowable cloud types in the original AS scheme, many clouds are sim- 
ulated frequently with each modifying the large scale environment some fraction of the total 
adjustment. The parameterization thereby “relaxes” the large scale environment towards 
neutrality. In addition to the RAS cumulus convection scheme, the GEOS-1 GCM employs 
a Kessler-type scheme for the re-evaporation of falling rain (Sud and Molod, 1988). The 
scheme accounts for the rainfall intensity, the drop size distribution, and the temperature, 
pressure and relative humidity of the surrounding air. 

Super-saturation or large-scale convection is defined in the GEOS-1 GCM whenever the 
specific humidity in any grid-box exceeds its super-saturation value. The large-scale pre- 
cipitation scheme rains at super-saturation, and re-evaporates during descent to partially 
saturate lower layers in a process that accounts for some simple micro-physics. 

The GEOS-1 GCM turbulence parameterization consists of elements which handle vertical 
diffusion (Helfand and Labraga, 1988) and surface fluxes of heat, moisture and momentum 
(Helfand, et al, 1991, and Helfand and Schubert, 1994). The vertical regime is divided 
into a free atmosphere, a surface layer, and a viscous sub-layer above the surface roughness 
elements. The turbulent eddy fluxes are calculated using a variety of methods depending 
on the vertical location in the atmosphere. 

Turbulent eddy fluxes of momentum, heat and moisture in the surface layer are calculated 
using stability-dependent bulk formulae based on Monin-Obukhov similarity functions. For 
an unstable surface layer, the chosen stability functions are the KEYPS function (Panofsky, 
1973) for momentum, and its generalization for heat and moisture. The function for heat 
and moisture assures non-vanishing heat and moisture fluxes as the wind speed approaches 
zero. For a stable surface layer, the stability functions are those of Clarke (1970), slightly 
modified for the momentum flux. The moisture flux also depends on a specified evapo- 
transpiration coefficient. 

Above the surface layer, turbulent fluxes of momentum, heat and moisture are calculated 
by the Level 2.5 Mellor-Yamada type closure scheme of Helfand and Labraga (1988), which 
predicts turbulent kinetic energy and determines the eddy transfer coefficients used for 
vertical diffusion. 
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4 The Analysis Scheme 


As mentioned above, the forcing terms for the IAU assimilation are obtained using a stan- 
dard statistical analysis methodology. A short term 6-hour forecast to the analysis time, 
t k (= 00, 06, 12, or 18 UTC) is used as the background, w£ (refer to Fig. 2)* Using all 
observations, w£, falling within a 6-hour data window centered on the analysis time, the 
analysis, w £, is defined as the background forecast plus a correction expressed as a quasi- 
linear combination of the differences between the observations and a, priori estimates of the 
observations computed from the background forecast, 

= wjj + £ fc (w£- U k w(). (1) 

In GEOS-1 DAS the “forward problem” transformation 'H k is factored as H k = H k S where 
S transforms from the model’s variables (temperature, wind, specific humidity, surface 
pressure) and sigma coordinates to the analysis variables (geopotential height, wind, water 
vapor mixing ratio, sea level pressure) and pressure coordinates. H k is the multi-linear 
interpolation of the pressure-level background vector w{ = £(w£) to our best a priori 
estimate of w£. Consistent with the factorization of the weights IC k also factor, fC k = 
$+K k , where K k is the usual gain matrix (see eq. 5) used to define the increments for the 
analysis variables on the constant pressure analysis levels, and S+ is the transform of these 
increments back to the model’s variables and sigma coordinates. The notation here follows 
that of Cohn and Parrish (1991). 

The forcing vector £ k used to modify the model’s tendencies during the IAU assimilation 
from time t k — 3 to time t k + 3 is held constant during that 6-hour interval and is obtained 
by transforming the pressure-level analysis increment back to the model’s sigma coordinates 
and variables (see Takacs et al. 1994 for details of the vertical interpolation). Thus & is 
defined by 

6 = (2) 

where N is the number of model time-steps of length A t seconds in the 6— hour assimilation 
update period centered on time tk. As mentioned above, S + represents the pressure to 
sigma transform for analysis increments. Since there are more model sigma levels than 
analysis pressure levels, the transform S is not actually invertible. Hence, S+ only denotes 
an approximate inverse for S. It does, however, act in such a way that = 0 in the interior 
of any region where Kfc(w£ — = 0, i.e., where the analysis increment is zero because 

there are no nearby data, the relevant components of will also be zero. 

By minimizing the variance of the analysis error, e£ = <S(w{, — w£), we are led to a system 
of linear equations which depends on the expected error characteristics of the background 
and observations and has the gain matrix Kk as its solution (cf, Jazwinski 1970, Example 

7.5). 

K k = [H Hi + -h] [H t H "I + «* + (H t J t + Jl Hi)] ■ ( 3 ) 

Using Kk from (3) in (1), we also obtain an equation for the analysis error variance, 

diag{f£} = diag{(7 - K k H k )P{ - KJl). (4) 
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In (3) and (4), Pg = (P k l a k T ) = SP£S T is the pressure-coordinate analysis error covariance 
matrix, P/ = T ^ = SPj[S T is the pressure-coordinate background error covariance 

matrix, Rk = is the observation error covariance matrix, and J* = is the 

pressure-coordinate cross-covariance matrix of background and observation errors. The 
background error is defined by e £ = *S(w£ — w£), and the observation error is defined by 
€ k = w 2 “ where w[ = <Sw£ is the unknown true state projected onto the analysis 

grid. We have made the assumption that both the background and observation errors are 
unbiased, i.e., that ^ i = (c£) = 0. The standard assumption will also be made that e £ 

and are uncorrelated, so that Jk — 0 (see Daley 1991, for the more general case). Given 
this last assumption (3) and (4) simplify to 


K k = [H k PlHl + R k ]~\ 

( 5 ) 

diag{F£} = diag{(/ - K k H k )P*}. 

( 6 ) 


The computational flow of the assimilation is as follows: 


Step 1. The innovation vector, [w£ - is computed, and the quality control of the 

observations is performed. See section 5 for details. 

Step 2. The components of w£ are partitioned into small geographically localized subsets 
of related variables. For example, one of these subsets might contain the height and 
wind components for the analysis grid-points in a small latitude-longitude-pressure 
box. For each of these subsets of w£ components, a local data selection is made from 
the quality checked innovation vector. See section 6 for details of this partitioning 
and the local data selection algorithm. 

Step 3. The innovation covariance matrix, [PfcP/ZfJ + #*], is approximated locally using 
the types and locations of the selected data. Equation (5) is then solved for those 
observation weights in Kk that give the contributions of the selected data to the w£ 
components being computed. Sections 7 and 8 contain details of the approximations 
for Rk and respectively. 

Step 4. The weights from Kk computed in step 3 are also used in (6) to estimate the 
analysis error variance components in diag{P£}. 

Step 5, The forcing vector £* is computed as in (2) and passed to the model for use in the 
assimilating IAU integration, which starts at time £* — 3 hours and extends to time 
tk + 3 hours. At tk +3 hours, £* is discarded and the integration is continued without 
analysis forcing to to produce the background, w £ +1 , for use in the next analysis. 

The forecast error variances, diag{P^ +1 ), which are also needed for the next analysis, 
are estimated using an error growth model together with the analysis error variance 
estimates for time tk that were computed in step 4. See section 8 for details of the 
forecast error growth model. 


These steps are further described in the following sections. 
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5 


Quality Control of Observations 


5.1 Pre-analysis Data Checks 

The GEOS-1 DAS ingests the global conventional observations and the temperature re- 
trievals from the HIRS2/MSU/SSU sounders on the NOAA satellites. Either the NESDIS 
retrieved temperature profiles or those created by the GLA physical retrieval system are 
used. Conventional data in NMC’s Office Note 29 , format as well as the NESDIS format 
retrieval data are unpacked and put into common format data sets of one day each. A day 
corresponds to the four analysis times: 0, 6, 12 and 18Z. Thus, the observations in one file 
will be from 2100Z of the previous day to 2059Z of the current day. 

The unpacking process keeps all data that appear in the original data sets, except those 
that do not have a realistic time stamp. Reports that are obviously in the wrong synoptic 
time (late arriving data, for example) are moved to the correct file. These data sets form a 
complete set of historical observations that are easily manageable on the computer system. 

A second preprocessing step standardizes the observations for ingest into the 01. Only 
observations of quantities to be analyzed are extracted for these data sets: sea level pressure 
and wind, upper-air height, wind and moisture. The satellite temperature retrievals are 
converted into thicknesses. The observations are stratified by type (e.g., surface land, 
surface ship) and in some instances by location. This is necessary to assign the proper 
observation errors to the data. More stringent checks on the data are performed: 

1. Quality marks provided with the observations are used to eliminate bad data outright. 
Data marked as being of suspicious quality are kept along with the quality flag. 

2. Observations with grossly bad values are deleted. 

3. A hydrostatic check is performed on rawinsonde data. 

4. Satellite profiles are checked for completeness. 


Before the observation data sets are ingested by the analysis scheme, a program is run to 
detect gaps in the observation time series. This program categorizes each observation by 
type, synoptic time and in the case of retrieved satellite temperature profiles, by the vertical 
extent of the profile ( stratospheric and tropospheric profiles in the case of NESDIS TOVS 
data). These counts are then graphically displayed as a series of bar charts as in Fig. 3. In 
this figure, the number of NESDIS TOVS path A (clear column) temperature retrievals for 
each synoptic time in the month of August 1985 is shown. The program has identified a 
significant gap in the TOVS data record from 11 to 17 August. In this case, it is likely that 
the data were not received properly from the primary archive center. This particular gap 
was not filled for the reanalysis, but usually such gaps have been filled with additional data 
received through secondary sources. These charts also serve as a record of the observational 
data that are available to the analysis at a particular time and are made available to the 
users. 
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Figure 3: Number of NESDIS TOVS path A (clear column) temperature retrievals for each 
synoptic time in the month of August 1985. Stratospheric retrieval profiles are those that 
only have levels above 100 hPa; tropospheric retrieval profiles are those that only have levels 
below 100 hPa. Complete profiles report at all levels. 
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5.2 Analysis Data Checks 


The data quality control is an important part of any data analysis scheme, as demonstrated 
in Shaw et al. (1986). More recent studies performed with the GEOS-1 DAS confirm 
that changes in quality control and subsequent data selection can have a major impact on 
the assimilation. The quality control technique currently employed in the GEOS-1 DAS 
consists of two major steps: a gross error check and a buddy check (Seablom, 1990). The 
gross error check is defined to be: 

A 2 < ((a 0 ) 2 + (a') 2 ) r, (7) 

where A is the difference between an observation and the interpolated background first- 
guess value, (a 0 ) 2 and {a f ) 2 are the observation and forecast error variances, respectively, 
and T is a subjectively defined tolerance value which varies with quantity, latitude and 
height. 

The tolerance value is somewhat reduced for all quantities in the tropics and increased 
slightly for the winds near jet level in the middle latitudes. Those data that fail to satisfy 
(7) are marked as suspect The buddy check involves performing a single pass successive- 
correction analysis of the data that passed the gross check to the locations of the suspect 
data. The difference between the interpolated value and the suspect value is then compared 
to the error statistics as in (7) and a decision is then made to reaccept the observation or 
to reject it. Typical rejection rates are between 5 and 10%. 


6 Localization of the Analysis Problem 


The GEOS-1 DAS global analyses are performed as a series of localized analyses on smaller 
regions referred to as mini-volumes . These mini-volumes are a set of non-overlapping groups 
of analysis grid-points. Associated with each mini-volume is an approximately cylindrical 
search region , having 3200 km diameter, from which the data are selected for the generation 
of the analysis in the mini-volume. The GEOS-1 DAS analysis has three distinct types of 
mini- volumes, each containing a different number of horizontal grid-points. In the trop- 
ics and low latitudes, each mini-volume contains six horizontal grid-points, rectangularly 
arranged. The mid-latitude region, between 30 degrees and 82 degrees latitude, contains 
mini-volumes with eight horizontal grid-points, while the polar regions place an entire lat- 
itude band of grid-points into each mini-volume. Additionally, each mini-volume contains 
two vertical layers of grid-points. The total analysis contains nearly 12,000 mini-volumes. 

A data search within each search region selects the 75 observations closest to the volume 
midpoint, with at most 60% of these from rawinsondes. The covariance matrix is then 
formed and the linear system of equations in (5) is solved for the weights of the observations 
Kk< With the current configuration of mini- volumes, neighboring search regions have an 
approximate 85% overlap. Because only 75 observations are selected, however, the actual 
overlap could be substantially less. 

The mini-volume approach to performing the global analysis is preferable to a single grid- 
point approach for two reasons. First, solving a local problem centered on a region rather 
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than on a particular grid-point significantly reduces the number of covariance matrices that 
need be set up and solved. The data search, also computationally expensive, is minimized 
by avoiding redundant calculations. Second, the independent nature of the mini-volumes 
allows the GEOS-1 DAS analyses to operate well in a parallel computing environment. The 
major drawback is that the limitation of 75 observations per mini-volume reduces the search 
region overlap, particularly in areas of high data density, leading to the adverse effect of 
boxiness in the analysis increments. 


7 Specification of Observation Error Covariances 


The elements of Rk = R = [r t J ] are the observation-observation error covariances, t.e. 
r tJ = (^e° x e° T ^ = r Jit . The time index A: will now be omitted, and subscripts will be 
used to denote vector and matrix elements. We write r t}J = where <rf 2 = 

a j 2 = ( e J 2 ) 1 anc * observation-observation error correlation. In the GEOS-1 DAS 

implementation a° and are obtained from tabulated values which vary with instrument 
type and level. Vertical interpolation using the tabulated values is done as needed. The 
observation error standard deviations for the instrument types being used in the DAO’s 
baseline five-year assimilation are shown in Table 1. 

The error standard deviations shown in Table 1 are in general larger than the standard errors 
in the measurements. The tabulated values contain an error of representativeness in addition 
to the instrument error. If w x (true) is the theoretical result of a perfect measurement 
corresponding to the real measurement w°, then [wj* — v/°(true)] is the measurement or 
instrument error. For the analysis, we have defined the corresponding observation error by 

e° = w? - (tfw*), = [w? - w° (trite)] + [w, ’(trite) - (i/w*),j . (8) 

The last term in brackets in (8) is what is called the error of representativeness (Lorenc 
1986). It depends on such things as grid resolution, subgrid-scale variability, footprint 
size for satellite instruments, and importantly on errors in the formulation of the “forward 
problem” operator H . These effects are not explicitly modeled in our system. The values in 
Table 1 have simply been inflated to account in a crude way for the representativeness error 
term. It should also be noted that observation error correlations for satellite soundings are 
probably in large part a reflection of the correlation between the representativeness errors 
for the observations. 

In general, the correlation, /i® , between e° y and e° is set to <5 tJ , the Kronecker delta, so that 
errors from distinct observations are assumed uncorrelated. However, in two situations 
is modeled as the product of an isotropic horizontal correlation, which depends 

only on the horizontal distance, s tJ , separating the observations, and a vertical correlation, 
v °' v , which depends only on their vertical separation. The first case occurs when the data 
are two height observations from the same radiosonde ascent. Here /i° }h (s v ) = 1 and v 0yV 
is similar to the function from Lonnberg and Hollingsworth (1986). Our current i/° fV has 
been extended in its vertical extent and its eigenstructure has been slightly adjusted, as 
shown in Table 2. The second case occurs when the data are two height observations from 
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Table 2: The vertical correlation of observation errors for two geopotential height observations from 
the same radiosonde ascent. 


Pressure level (hPa) 


level 20 30 


20 

1.00 

0.00 

30 

0.00 

1.00 

50 

0.00 

0.00 

70 

0.00 

0.00 

100 

0.00 

0.00 

150 

0.00 

0.00 

200 

0.00 

0.00 

250 

0.00 

0.00 

300 

0.00 

0.00 

400 

0.00 

0.00 

500 

0.00 

0.00 

700 

0.00 

0.00 

850 

0.00 

0.00 

1000 

0.00 

0.00 


50 

70 

100 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

1.00 

0.89 

0.88 

0.89 

1.00 

0.89 

0.88 

0.89 

1.00 

0.84 

0.86 

0.88 

0.81 

0.83 

0.86 

0.73 

0.75 

0.78 

0.70 

0.71 

0.73 

0.53 

0.56 

0.59 

0.43 

0.46 

0.51 

0.24 

0.30 

0.36 

0.01 

0.04 

0.10 

0.00 

0.01 

0.03 


150 

200 

250 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.84 

0.81 

0.73 

0.86 

0.83 

0.75 

0.88 

0.86 

0.78 

1.00 

0.88 

0.81 

0.88 

1.00 

0.84 

0.81 

0.84 

1.00 

0.76 

0.76 

0.84 

0.62 

0.61 

0.71 

0.54 

0.53 

0.57 

0.45 

0.46 

0.40 

0.23 

0.24 

0.14 

0.07 

0.08 

0.03 


300 

400 

500 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.70 

0.53 

0.43 

0.71 

0.56 

0.46 

0.73 

0.59 

0.51 

0.76 

0.62 

0.54 

0.76 

0.61 

0.53 

0.84 

0.71 

0.57 

1.00 

0.78 

0.61 

0.78 

1.00 

0.78 

0.61 

0.78 

1.00 

0.39 

0.54 

0.70 

0.13 

0.14 

0.24 

0.05 

0.09 

0.15 


700 

850 

1000 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.24 

0.01 

0.00 

0.30 

0.04 

0.01 

0.36 

0.10 

0.03 

0.45 

0.23 

0.07 

0.46 

0.24 

0.08 

0.40 

0.14 

0.03 

0.39 

0.13 

0.05 

0.54 

0.14 

0.09 

0.70 

0.24 

0.15 

1.00 

0.64 

0.36 

0.64 

1.00 

0.48 

0.36 

0.48 

1.00 


the same TOVS sounding or from different soundings of the same type. In this case, the 
vertical correlation (see Table 3) is also a slight modification of that used in B87. The 
horizontal observation error correlation between NESDIS TOVS soundings of the same 
type is modeled by p 0 ’ h (s tJ ) = exp [— (s,j/150.0) 2 ], where s tJ is the separation in km. This 
function is considerably sharper than the functions which appeared in B87. If the soundings 
are of different types, then /if = S tJ is used. 


8 Specification of Background Error Covariances 

The central feature of Kalman filter theory (Jazwinski 1970, Daley 1991, Cohn 1993, and 
references therein) is its inclusion of a prediction equation for forecast error covariances, 

= +Qk-l, (9) 

where A k - 1 represents a linearization of the atmospheric prediction model and Q k -i is the 
covariance matrix of model errors, i.e., Q k _ x with wj. = + e t k _ l . 

Optimum interpolation can be viewed as an approximation to the Kalman filter in which 
(9) is replaced by an empirical model for the evolution of P/ (see Todling and Cohn 1994). 

Note that (9) gives some guidance in the design of such empirical models in that it shows 
that P/ depends upon: (1) the recent history (contained in F£_,) of the distribution and 
quality of assimilated observations; (2) the model error covariance (contained in Q k ~ i); 
and (3) the dynamical evolution of existing errors (governed by .4*_i). In the GEOS-1 
DAS implementation of 01, the last of these three factors is ignored and the first two are 
accounted for only in a crude way. Writing the elements of P/ = | J as p { } = we 

separately model the error standard deviations, <j{ and <jj , for the tth and jth components 
of the background error vector, and the correlation, between these error components. 
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Table 3: Vertical correlations of observation errors for two geopotential height observations from 
NESDIS TOVS retrievals of the same type. The correlation function from the upper triangular part 
of the table is used for type A (clear) and type B (partially cloudy) retrievals. The function from 
the lower triangular part of the table is used for the type C (cloudy) microwave retrievals. 


level 






Pressure level (hPa) 






20 

30 

50 

70 

100 

150 

200 

250 

300 

400 

500 

700 

850 

1000 

20 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

30 

0.00 

1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

50 

0.00 

0.00 

1.00 

0.87 

0.82 

0.67 

0.57 

0.45 

0.31 

0.28 

0.26 

0.22 

0.12 

0.00 

70 

0.00 

0.00 

0.87 

1.00 

0.83 

0.70 

0.58 

0.45 

0.30 

0.26 

0.23 

0.17 

0.06 

0.00 

100 

0.00 

0.00 

0.81 

0.85 

1.00 

0.83 

0.70 

0.57 

0.43 

0.38 

0.35 

0.25 

0.16 

0.00 

150 

0.00 

0.00 

0.68 

0.73 

0.83 

1.00 

0.81 

0.71 

0.59 

0.53 

0.48 

0.33 

0.22 

0.00 

200 

0.00 

0.00 

0.56 

0.61 

0.73 

0.84 

1.00 

0.85 

0.77 

0.65 

0.55 

0.28 

0.15 

0.00 

250 

0.00 

0.00 

0.46 

0.50 

0.63 

0.75 

0.85 

1.00 

0.86 

0.77 

0.65 

0.32 

0.16 

0.00 

300 

0.00 

0.00 

0.35 

0.39 

0.53 

0.65 

0.79 

0.86 

1.00 

0.84 

0.73 

0.37 

0.19 

0.00 

400 

0.00 

0.00 

0.34 

0.36 

0.47 

0.55 

0.67 

0.77 

0.84 

1.00 

0.85 

0.56 

0.36 

0.00 

500 

0.00 

0.00 

0.29 

0.31 

0.50 

0.47 

0.55 

0.64 

0.72 

0.85 

1.00 

0.72 

0.52 

0.00 

700 

0.00 

0.00 

0.24 

0.24 

0.30 

0.31 

0.32 

0.38 

0.47 

0.67 

0.81 

1.00 

0.82 

0.00 

850 

0.00 

0.00 

0.16 

0.12 

0.19 

0.20 

0.17 

0.20 

0.29 

0.51 

0.67 

0.83 

1.00 

0.00 

1000 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

1.00 


8.1 Background Error Variances 


The background error standard deviations a{ for geopotential height, sea level pressure, 
and water vapor mixing ratio are estimated in a two-step process. First, the analysis error 
standard deviations, <r, a , for these quantities at time tjt_i (square roots of elements along 
the diagonal of Pg_ j) are computed using a localized version of (6). A growth term is then 
added to simulate the effect of Qk-i °n the diagonal of Pj. (Pfaendtner and Sivakumaran, 
1991). The equation used for o{ at time t* can be written 


a{ = < 


+ 


growth term = + <*{l — 


( 10 ) 


The parameters a and 0 used in the definition of the growth term, which determine the 
growth rate and saturation values for the errors, respectively, depend on the pressure level. 
For geopotential height and sea level pressure there is an additional crude dependence on 
latitude. This latitudinal dependence is defined using tropical, £*< r and 0t r , and extra- 
tropical, a x t and Art. parameter values with a smooth transition between 15° and 45° north 
and south latitude. For example, for 0 we have: 


0 = 0{<p) = I 


0tr i 

\{0tr + At) + |(Ar — Art) cos (^25 
Art i 


M < 15° 

15° < M < 45° (11) 

45° > M . 


The form of this functional dependence on latitude was motivated by an examination of the 
growth of zonally averaged RMS errors from an ensemble of short-term (to 48 hr.) forecasts. 
The parameters (see Table 4) have been empirically adjusted by comparing the mean square 

innovation vector components, / |w" - (iTwA),J Y to the 01 estimate (cr °) 2 + ^( Ha ^),^ . 
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An example of the effectiveness of the error growth model is shown in Fig. 4. The figure 
shows a time series of 500 hPa of forecast error standard deviation estimated from the 
innovations and observation error standard deviations by the formula 

a s w((w°- Hw f )i - (O 2 ) 172 . (12) 

compared to the forecast error standard deviation produced using the growth model, eq. 
(10). This figure shows modelled and observed a 1 for a box in the South Pacific (165W- 
105W ,60S-36S) , from July 7 to August 4, 1979. Although the magnitude of the modelled a 5 
is slightly larger than observed, both curves are highly correlated after July 14. Therefore, 
the empirical growth model seems to capture the main characteristics of the temporal 
variations of the forecast error standard deviations. 



Date 


Figure 4: Upper panel: Observed (solid line) and modelled (dashed line) 500 hPa forecast error 
standard deviations for a box in the South Pacific (165W-105W, 60S-36S) from July 7 to August 4, 
1979. Lower panel: the number of observations used to estimate the observed forecast error standard 
deviations. 


Outside of the tropics the standard deviations of the u and v components of the wind errors, 
o™ and cr x , used in the multivariate upper air height-wind analysis are determined using 
the geostrophic approximation as in B87 (see also Appendix A), 


K ) 2 = «) 2 = 


g 2 «) 2 f d 2 n z ' z {s) \ •> 
|/| 2 l 8s 2 UoJ’ 


(13) 


where g is the gravitational constant, / the Coriolis parameter, a z x the forecast height error 
standard deviation from (10), and /i z,z (s) the isotropic separation dependent model function 


15 




Table 4: Parameters used in the background error variance growth model for geopotential height, 
water vapor mixing ratio and sea level pressure. 


level 


Geopotential height (m) 


Mixing ratio (g/kg) 

Tropics 

a tr 

Extra-Trop. 

Oxt 

Tropics Extra-Trop. 
Ar Pxt 

a 

0 

10 

125.00 

70.00 

147.00 

112.50 

— 

— 

30 

110.00 

59.00 

132.00 

109.50 

— 

— 

50 

95.00 

48.00 

117.00 

106.50 

— 

— 

70 

80.00 

37.00 

105.00 

103.00 

— 

— 

100 

60.00 

35.00 

96.00 

98.50 

— 

— 

150 

25.00 

32.00 

76.00 

95.00 

— 

— 

200 

21.00 

30.00 

56.00 

91.50 

— 

— 

250 

17.00 

28.00 

37.00 

88.50 

— 

— 

300 

14.00 

26.50 

29.00 

85.00 

0.05 

0.40 

400 

13.00 

23.00 

26.90 

78.00 

0.15 

1.20 

500 

12.00 

20.00 

24.50 

71.00 

0.20 

1.60 

700 

10.50 

15.00 

20.00 

58.00 

0.65 

5.20 

850 

9.00 

12.50 

20.00 

56.50 

0.90 

7.20 

1000 

9.00 

10.50 

20.00 

55.00 

1.20 

9.60 

SLP 

1.13 

1.31 

2.50 

6.88 

— 

— 


for the horizontal correlation of forecast height errors, which is discussed in the next section. 
In the tropical band 15°S < f < 15°N, a" and are assigned the pressure level dependent 
value given in Table 5. At points in the subtropical bands, 25°S < <p < 15 S and 15 N 
< <p < 25°N, a“ and a ? are defined by 


K ) 2 = K ) 2 = 


f25-Ml 

(rr u ' v l 2 4- 

M - 15' 

s 2 «) 2 

d 2 fi z,z (s) 


h- * 

O 

V^trop' ' 

10 J 

I/ 25 I 2 

ds 2 

0- 


(14) 


where / 2 s = 2flsin(25°). 


Outside of the tropical band, 20°S < cp < 20°N, the standard deviations of the u„ and v, 
components of the near surface wind errors over the oceans are defined as in B87 by 


K-) ! = W ) 2 = ^(C? + C, ! ) 


d 2 /z p - p (s) 


ds 2 

5 = 0 - 


(15) 


where is the background sea level pressure error standard deviation from (10), a the 
Earth’s radius, C, and F, the coefficients for the surface wind model derived in B87, and 
/z p,p (s) the isotropic separation dependent model function for the horizontal correlation of 
background sea level pressure errors discussed in the next section. In the tropical band, 
10°S < <p < 10°N, cr“ J and cr, Va are set to 5 ms -1 , while in the subtropical bands, 20°N 
< ip < 10°S and 10°N < <p < 20°N, they are defined by 


«*) 2 = {<•? = 




( 16 ) 
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Table 5: Tropical upper air wind forecast error standard deviations. 


Tropical 

Pressure interval 

<t u and <x v (ms" 1 ) 


P < 

20 

7.1 

20 

< P < 

40 

7.1 

40 

< P < 

60 

8.0 

60 

<P< 

85 

7.3 

85 

< P < 

125 

7.3 

125 

< P < 

175 

7.3 

175 

< P < 

225 

7.4 

225 

< P < 

275 

7.0 

275 

< P < 

350 

6.2 

350 

< P < 

450 

5.2 

450 

< P < 

600 

4.5 

600 

< P < 

775 

4.4 

775 

<P< 

925 

4.4 

925 

< P 


4.3 


8.2 Background Error Correlations 


For all permutation pairs, [a, 6], of variables for the upper-air multivariate analysis [geopo- 
tential height ( 2 ), east-west wind component (u), and north-south wind component (t?)] the 
forecast error correlations, are assumed separable i.e. they are written as the product 
of a horizontal, p a ' b (s tJ ), and a vertical, t'®- 6 , error correlation. As in B87, the horizontal 
height-height error correlations are modeled using a damped cosine function (Thiebaux 
1975), 

/**’*(«) = - X , [ci cos(c 2 s) -I- c 3 ][l + (c 4 s) 2 ]- c *. (17) 

Notice that this function has two length scales, c 2 and c 4 . The fitting parameters have been 
recomputed since B87, and a single function is used globally in contrast to the regionally 
dependent functions described in B87. The current parameter values (cj = 0.0129928, 
c 2 = 0.00389265 km" 1 , c 3 = 0.694005, c 4 = 0.00105123 km" 1 and c 5 = 1.20815) were 
obtained by a fit using 500 hPa radiosonde height data from North America for January 
and February of 1979. The graph of (cj + c 3 )p z ’ z (s), illustrating its fit to the data, is 
shown in Fig. 5. The vertical correlation of background height errors, v z f , is obtained by 
interpolating from the values appearing in the upper triangular part of Table 6. In GEOS-1, 
the chord length approximation is used to computed distances (see Appendix 9). 


For the correlation of forecast height errors with forecast wind errors as well as the wind- 
wind forecast error correlations an approximation to the geostrophic assumption is used. 
Starting with the height-height error covariance model, = <r z a z fi z ’ z v z j z , and the 

geostrophic approximations for the wind component errors, 


e u = - 


g_d£ 
af d<p 


and 


9 de z 
af cos <p dX ’ 
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(18) 



Table 6: Vertical correlation of background errors for geopotential height (shown in upper triangular 
part of table), and for water vapor mixing ratio (shown in lower triangular part of table). 


level 






Pressure level (h Pa ) 





1000 

10 

30 

50 

70 

100 

150 

200 

250 

300 

400 

500 

700 

850 


1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

30 


1.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

50 



1.00 

0.87 

0.71 

0.29 

0.19 

0.06 

0.00 

0.06 

0.08 

0.10 

0.14 

0.09 

70 




1.00 

0.75 

0.31 

0.19 

0.05 

0.00 

0.05 

0.06 

0.09 

0.12 

0.09 

100 





1.00 

0.60 

0.46 

0.30 

0.21 

0.22 

0.21 

0.15 

0.18 

0.12 

150 






1.00 

0.82 

0.66 

0.56 

0.47 

0.43 

0.25 

0.19 

0.10 

200 







1.00 

0.81 

0.72 

0.61 

0.52 

0.32 

0.24 

0.14 

200 








1.00 

0.86 

0.77 

0.71 

0.49 

0.40 

0.25 

300 









1.00 

0.85 

0.80 

0.59 

0.48 

0.30 

400 









0.94 

1.00 

0.87 

0.70 

0.59 

0.37 

500 









0.60 

0,71 

1.00 

0.77 

0.68 

0.43 

700 









0.07 

0.12 

0.33 

1.00 

0.83 

0.57 

850 









0.02 

0.04 

0.13 

0.50 

1.00 

0.66 

1000 









0.04 

0.05 

0.04 

- 0.01 

0.27 

1.00 


a fj 


d V ' 


af, \V dtp. 


afj cos ifj 


g 

afj d>fj 


d\)) 


9 /(d? 

a f % cos \ \ dX 


)/*) <*/. 


we get: 

(^> 

(^) 

<^> 

(C?) 

(<7<7 ) 

( e v e v ) _ 

\ ’ - , / a 2 f t fj cos cosipj 

In computing the partial derivatives of /cf e j) in (19)-(25) we make the additional assump- 


g 2 (—) \ = g2 

a 2 /./, \ V / , \ dtp ) jf a 2 f t fj d<p t d<fij 


d<p t 



9 

d 


afj cos <Pj 

dXj 

9 

d( 

f e z e z ) 

3/ 

a/» cos <p t 


dX , 

9 2 

d 2 



((£ 




a 2 («'?*;) 


a 2 fJ 3 cos<pj\\dy> ) x \d\ ) J a 2 f t fjCos<pj d<p % d\ 3 


a 2 f,fjCos<p, 

„2 


cosy), d\ t d(pj 


(19) 

( 20 ) 
( 21 ) 
( 22 ) 

(23) 

(24) 

(25) 


((».(£)}-«* 

//a?', /a?\ \ = g 2 a 2 (?q) (26 . 

\ v dX ) , V dX ) jj a 2 f t fj cos y?, cos <p : dX x dX } 
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Figure 5: Fit of scaled damped cosine correlation function to radiosonde minus background 
500 hPa height difference correlations using North American radiosondes for January and 
February 1979. 


tion of local homogeneity for the height error variances, 


da z ^ do z 


0 . 


(27) 


Consequences of this homogeneity assumption are discussed in Cohn and Morone (1984) and 
Morone and Cohn (1985). Derivations of effective computational equations for the height- 
wind and wind-wind error correlations including modifications for the tropics are given in 
Appendix C. Note that separate methods have been used to decouple tropical wind-wind 
variances (eq. 14) and tropical wind-mass correlations (eqs. 118,120-123). 


9 Current Development Status of GEOS DAS 


In the time since GEOS-1 DAS was created to perform the 5 year reanalysis, one new 
version of the analysis component has been made operational (01 1.5) and another version, 
with many more significant changes, is about to become operational (OI 2.0). The most 
important change in version OI 1.5 is the correction applied to the wind- wind forecast 
error correlations in the deep tropics. The geostrophicaliy derived correlations are not used 
across the equator because sign changes introduce spurious divergence in the wind analysis; 
no changes were introduced in the extratropics (defined as those latitudes poleward of 15°). 

For OI 2.0, the code has been generalized to allow up to 300 observations (one obser- 
vation equals one observed variable) to contribute to the analysis of the grid points in 
each mini-volume. We have also added the ability to assimilate the SSM/I surface wind 
speeds, SSM/I total precipitable water retrievals, and have experimented with assimilation 
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of the microwave limb sounder temperatures from the Upper Atmosphere Research Satel- 
lite (UARS/MLS). The code has been generalized to accept new conventional data types, 
such as the ARINC Communication Addressing and Reporting System (ACARS) aircraft 
data, which are winds reported by wide-bodied US commercial aircraft. The large mini- 
volumes at the poles, which had been a collection of all pole points, have been broken up 
into mid-latitude sized mini-volumes. This eliminated some discontinuities associated with 
such large mini-volumes. The observation increment history dataset now contains a history 
of the quality control decisions, and we are carrying the time of each observation through 
the analysis to begin experimenting with the decoupling of satellite temperature retrievals 
from different orbits. Work is in progress to retune observation and forecast error statistics 
consistent with recent improvements in the system. 

Since GEOS-1 DAS, there have also been several improvements made to the GEOS GCM. 
The first improvement entails the development of a fourth-order version of the Sadourney 
scheme for the momentum equation. The “fourth-order Sadourney” is derived in Suarez 
and Takacs (1995). This scheme conserves total energy and potential enstrophy for non- 
divergent flow, and is fourth-order in the sense that the advection of the second-order 
vorticity reduces to that governed by the fourth-order Arakawa (1966) Jacobian. This 
scheme was chosen rather than the Arakawa-Lamb (1981) scheme modified to fourth-order 
by Takano- Wu rtele (1982) due to severe polar noise problems generated in three-dimensional 
GCM calculations from the terms required for potential enstrophy conservation for general 
flow. Together with the “fourth-order Sadourney” scheme for the momentum equations, 
the current GEOS GCM also uses fourth-order horizontal advection in the thermodynamic 
and moisture equations. It uses a scheme developed by Arakawa which has been used for 
many years in the UCLA GCM. 

The dynamical climatology of the GEOS GCM has been thoroughly examined by Takacs 
and Suarez (1995). They have shown that increasing the resolution and/or the order-of- 
accuracy in the GEOS GCM has a very significant impact on the zonal mean flow. Running 
the second-order GEOS-1 GCM at both 4° x 5° and 2° x 2.5° horizontal resolution, an 
examination of the first and second moments from a five year simulation revealed a strong 
systematic bias for most fields between the two resolutions, particularly in the southern 
hemisphere. Most of these biases are related to the simulation of the transient eddies, with 
increased accuracy enhancing transient transports of heat and momentum and decreasing 
transports of moisture. These systematic biases are removed, however, when fourth-order 
accuracy is used. The higher-order and higher- resolution experiments are in closer agree- 
ment to the GEOS-1 DAS analysis in the southern hemisphere where the transports are 
primarily determined by the transient flow. 

Another major enhancement to the GEOS GCM is the ability to run both simulations 
and the GEOS DAS assimilation using a rotated coordinate. In the current system, the 
geographic placement of the computational pole is arbitrary. We have found that due to 
variations of the scheme near the poles required from conservation constraints, the compu- 
tational instability discussed by Hollingsworth et al. (1983) for the nearest grid-point to the 
pole gives rise to polar noise when confronted with strong cross-polar flow. This instability 
is proportional to the mean zonal wind speed and the Coriolis parameter. By rotating the 
computational grid to the geographic equator, however, the instability near the computa- 
tional pole is removed due to the vanishing Coriolis term. In addition, the geographic pole 


20 


now using the transformed grid is also free of noise. 
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Appendix A: Chord-length/ Great Circle Distance Approxi- 
mation 


The background and observation error correlation models used in the analysis are expressed 
as isotropic functions of the separation distance between points on the sphere. Given two 
points Pi = (A,-, <pi) and Pj = (A j,ifij) on the sphere, we can write their position vectors as: 


p, = aq, = a [cos ipi cos A,, cos (p, sin A,-, sin <p,] 


(28) 


and 


Pj = aq j = a [cos <fj cos A j , cos <pj sin A j , sin <pj ] , (29) 

where a is the Earth’s radius. The true shortest distance arc length along the great circle 
connecting P, and Pj is given by 


s tJ = a arccos(q, • q_,). (30) 

On the other hand, the square of the length of the chord connecting Pj and Pj is: 

S ij = (Pi ~ Pj) ■ (p, - Pj) 

= 2a (1 — q« - qj). (31) 

In GEOS-1 DAS the chord length approximation (31) is used to compute distances. For 
the sake of completeness, this section will derive equations based on the great circle formula 
(30) and the chord-length formula (31). 


We will need the first and second partial derivative of Sij and S{j with respect to <p iy 
and Aj. Note that if 7 and 77 are any two (possibly the same) of these four independent 
variables we have 


^L S in(^) = ?pL (fii) = _ Q fl(q»-qj) 

dy K a dy ^ a dy 


(32) 


and 


Pfij 

dr)d y 


a sin(-^-) -f- cos(— ) 

CL & 


ds^j 

dr) dy 


In particular, we need the following partial 
error covariance computations. 


. ds i:) ds i} o d 2 {q t -qj) 

dydy 13 dr/ dy drjdy 


(33) 


derivatives for the height-wind and wind-wind 


• For ( zu ), 


• For ( uz ), 


• For ( zv ), 


k dsij - ds,j 

— = = ~ ac iy 

d<pj d?j y 


i @Sj j • ds,j 

b Wi = V = 


At = 

dXj dXj 


flC|T 
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(34) 


(35) 

(36) 


lMTgNllCNAl.LV BLANK 



• For (vz), 


^ 7 


: dsij 

b 8K = 

(37) 

• For (uu), 

.3 ^ 2s «j _ 

dwdcpj 

— 0,{b^ Cyy "I" C%jC{yCyj ) 

(38) 


- 3 d 2 S{j _ 

d<pid<pj 

a(b Cyy -f* CjyCyj) . 

(39) 

• For (uu), 

.3 _ 

a^5A, 

— U^b^Cyx + CijCyjCix ) 

(40) 


£3 5 2 S tJ _ 
d<pid x, 

— &(y Cy X "I - CyjCjjp). 

(41) 

• For (trn), 

3 5 2 S,j _ 

dXidifj 

— d{b 2 C X y 4- CijCxjCiy ) 

(42) 


i3 d 2 s*j 

dXid^pj 

^ A * 

— d(6 Cxy H" Cxj c ty)* 

(43) 

• For ( vv ), 

3 5 2 Sij _ 

dXidXj 

Cxi + CijCix^xj) 

(44) 


II 

<c»T ^ 

CM J* 

^ £ 

-a(b 2 c X x + CixCxj )* 

(45) 

Where we have used, 



b = 

s « n (^-) = ^i-4* 


(46) 

b = 

(^) = v/2(l - «i). 


(47) 

C «'j = 

q. • qj 


(48) 

— 

+ cos ipi cos A, cos ipj cos Xj + cos (fii sin A,- cos ifj sin Xj + sin </>, sin ipj, 

== 

dqj 

q ‘ ' dXj 


(49) 

= 

— cos (pi cos A, cos pj sin Xj + cos <pi sin A,* cos p>j cos Xj , 


= - cos <pi cos A; sin ipj cos A j - cos <pi sin A,- sin ipj sin Xj + sin y>,- cos <pj , (50) 

dq, 

Cxj ~ dx i ‘ qj 
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Appendix B: Damped Cosine Correlation Model 

As in Baker et al. (1987), the horizontal height-height error correlations, as well as the 
sea level pressure-sea level pressure error correlations, are modeled using a damped cosine 
function (Thiebaux 1975), 


n = n zz (s) = v pp (s) = — 1 — [ci cos(c 2 s) + c 3 ][l + C4S 2 ] C5 . (57) 

^1 r C3 

The fitting parameters have been recomputed since Baker et al. (1987), and a single function 
is used globally in contrast to the regionally dependent functions described in that paper. 
The current parameter values, 

ci = 0.0129928, (58) 

c 2 = 0.00389265 km -1 , (59) 

c 3 = 0.694005, (60) 

c 4 = 0.00105123 km -1 , (61) 

c 5 = 1.20815, (62) 


were obtained by N. S. Sivakumaran using a fit to 500 hPa radiosonde height data from 
North America for January and February of 1979. 

If we define 

(63) 

(64) 
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Ci = 


Cl 

C 1 + c 3 ’ 
C3 


C3 = 7 = 

C 1 + C 3 



11 i j.UIuii:i .infill LI ..iilil ll jl 


we can write 


p 2 = 

7 = 

*7 = 

A = 
A' = 
B - 
B' = 

C = 

C' = 

D = 


cic\ + 2C4C5, 

ci c\ 


P ’ 

2cjc 5 

P 2 


= 1- 7, 


(1 + c^s 2 )- 1 , 

-2c 2 A 2 s, 

A c >, 

- -T)0 2 ABi s, 
sin(c 2 s) 

(l/s)[cos(c 2 s) - C], 

1 „ 2 sA' 

2c * As = ~~A~~^B 


H = [cj cos(c 2 s) 4- c 3 ]jB. 


Differentiating we get 


p! — (fi/B)B' - ci c 2 sin (c 2 s)Z? 
= -r)f3 2 Aps - c\c\CBs 
= -P 2 s[rjAp + jBC], 


(65) 

( 66 ) 

(67) 

( 68 ) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 


(75) 


(76) 


and differentiating again gives 

p" = -p 2 (^ rs )-P 2 s[r 1 (A^ + A , f i) + y(BC , + B , C)] 

= -p 2 {v{Ap + sAp' + sA'n ] + 7 [BC + SBC' + sB'C ]} 

= -p 2 {AfJ,r)[l - 2c\As 2 {\ + c 5 )] + fl7[cos(c 2 s) - 4C4C5ACS 2 ]} 

= -/3 2 {ApT][l - (1 4- c 5 )D] + B7[cos(c 2 s) - 2c 5 CD]}. (77) 

At the origin we have: p(0) = 1, lim a _,o ( M ^ ) = ~P 2 i an ^ = ~P 2 ■ The graphs of 

the functions p, -V/ (sp 2 ), and -p"/P 2 are shown in Fig. 6. 


Appendix C: Height-Wind and Wind-Wind Background Er- 
ror Covariances with Damped Cosine Model Function 

In this appendix we give the explicit formulas for the correlation functions used in the 
multivariate wind-height analysis. 
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Damped Cosine and its Derivatives 



Figure 6: The functions p, —p' /(s/3 2 ), and —p"//3 2 calculated using the parameter values 
given in the text. 


C.l Correlation modeling assumptions for background errors 


For all permutation pairs, ab , of the three upper-air analysis variables [geopotential height 
( 2 ), east-west wind component (it), and north-south wind component ( v )] the correlation 
between background errors in variable a at location Pi = (A,-, (pi^pi) and background errors 
in variable b at location Pj = (Aj, <Pj,pj) is defined by 



° a i°) ' 


(78) 


Our first modeling assumption, often referred to as the separability assumption, is that each 
of the /*(e®£j)’s can be written as the product of a horizontal, pfj = t* ab (sij), and a vertical, 
v i j = > error correlation, where Sij = s(A,*, <£>,-, Aj, (pj) is a measure of the horizontal 

distance between points P t and Pj, and Sij — S(pi,pj ) is a measure their vertical separation, 

i.e., 


« M a6 (^>° 6 (%) = 


(79) 


The second modeling assumption is that the background wind component errors can be re- 
lated geostrophically to the background height errors, i.e. we will make use of the geostrophic 
approximations for the background wind errors, 


g 

o,fi d<pi 


and 


def 


a fi cos d\, 


(80) 


The geostrophic assumption implies that we can express all of the (c“ej)’s in terms of (efe?) 
and its partial derivatives with respect to A,-, <pi, A j and ipj: ' 




(TiCTjPijUij 


(81) 
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(82) 




<4# 

<«?$> 

<«T$> 


5 /f , £fl v _ g 


— (c* — ) = 

~ r \ c * . / 


d(efej) 


afj dvj 

9 / de i ,z\ g &(*!$) 

afid(fii * afi dfi 

9 - —l gMfH 

\ e * a\ . ' 


a /, cos ' 1 5Aj afj cos <pj dXj 

9 « _ g ^11 

J' 


a/t cos N $A t 


a/j cos 


g 2 ,de?dej = g 2 &(±J) 

a 2 fifj ' <V, 0y>j ' oPfifj d<pid<pj 


-( 


fc? *|v = £ 

A \ ^2 A T/ 




a 2 fifj cos (fj ' cV, 5Aj ' a 2 fifj cos ipj dtfidXj 

g 2 Mde = g 2 
a 2 fifj cosifi d\id<pj a 2 fifj cos y>; 3X,d<pj 


a? fifj cos (piCO&(pj^d\i d\j' a. 2 fifj cos <p,- cos <pj dXidXj 


> = 


gj^g 


(83) 

(84) 

(85) 

( 86 ) 

(87) 

(88) 

(89) 


The third modeling assumption is that the height error standard deviation is horizontally 
homogeneous. This homogeneity assumption (see Morone and Cohn 1985) is expressed 
through the approximations, 

3o do ron'i 

(90) 

and will be used to simplify the partial derivatives of which appeared in (82) through 

(89). 


C.2 Derivation of wind component error variances 

In order to derive expressions for a u and a v which are consistent with our correlation 
modeling assumptions, we will first examine ) and (eJ'Cj) in detail. From the separability 
assumption we have: 

( 91 ) 

Using our geostrophic approximation and separability a second time, this becomes 

~ g 2 

, ~ a 2 fjj 3<pid<pj 


g2 ^(gf^/if/lff) 

~ a 2 fifj 3<pidipj 

With two applications of the homogeneity assumption, we arrive at 


(92) 


a x a j Pa v a ~ 


g 2 d 


g 2 **f^/y* 
a 2 /.7j ‘ u 


(93) 
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( 94 ) 


We now use the chain rule to rewrite the expression in brackets to get 

tT u rr u n uu i/ uu ~ „z _z f ® Pij ds ds , d 2 S } zz 

' 3 ^' 3 ' 3 ~ a 2 fifj ’ J * ds 2 d<pid<pj + ds d<p,d<p : J ij ' 
Taking the limit as j -4 i, or equivalently as S — ► 0, we have 


^2 = g 2 (<rf) 2 f^V/ 5s as + 5/ifJ 5 2 s~ 




a 2 /, j-m 1 


5s 2 5y>, 5y?j 5s difidipj 


;}• 


(95) 


^ t*ij i s modeled with the damped cosine function described in the previous section, we can 
use equations (75), (76), (34) and (35) to simplify and get 


K ) 2 = - 


7 2 /? 2 K) 2 

2A2 


a 2 fi 


Hm r^w 

i-H'l fc2 /’ 


b 2 


(96) 


where /3 2 depends on the damped cosine fit parameters as given in (58)— (62); c iy and c w 
are defined in (50) and (52); and b = b from (46) if s is the true arc length, or b = 6 from 
(47) if s is the chord length distance approximation. In any case, the limit in (96) is equal 
to — a 2 and we have: 


Jx 


(97) 


as appeared in Baker et al. (1987). This clearly needs to be modified near the equa- 
tor. These modifications, the so-called decoupling of the geostrophic approximation in the 
tropics, will be discussed in subsection C.4 below. 


The analog of (94) for (cfc/) is 


r V -V , t vv , ,vv 




,/dV?; ds ds 5/if/ d 2 s 


+ 


a 2 fifj cos (pi cos ipj ' u ds 2 5A, d\j ' ds dX,dX : 
Taking the limit a s j i and using equations (75), (76), (36) and (37) we get: 


(98) 


W ) 2 = \m{ T a ' c “ Cli }, 

a fi 1 b 2 cos cos<£, > 


COS (fi cos <fij 


(99) 


where c, T and c x j are defined in (49) and (50). The limit in (99) is also equal to —a 2 . 
Hence, (erf) is the same as (erf) 2 as given in (97). Notice that <r“ equals a v only under the 
assumption of homogeneous variances expressed by (90). 


C.3 Formulas for computing background error correlations for the upper- 
air analysis 

Now that we have formulas for the background wind error variance, we can substitute 
them into our separable background error covariance models to get effective formulas for 
computing the various error cross-correlation functions. We start with /if/*. 
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Substituting erf and erf as given by (97) into (94), and noting that vf? must equal v- , we 
get: 

sign(y,) sign(y?j) ( d 2 tfj ds ds djJ^_ \ ( 100 ) 

n 2 f3 2 l ds 2 dip, d<Pj ds d<fiid<Pj ' 

Substitution using (73) and (74) allows us to rewrite this as 

, lf s rMi*'( 5 )i d 2 s 1 [**“"(«) 1 \ nrm 

/ iV; = sign(^sign(^){-^[^^]^^--[ _^ 2 \CiyCyj)- (101) 

If we are using the chord length approximation for horizontal separation, i.e., s s, we can 
substitute using (39) to obtain, 




= ngn(v»i)8ign(v>i)|[-^ja| 


Cyy 4" 


1 


^ ^ (S/a) 2 

= sign(vJi) sign(<pj){.R(s)c yy + T(s)c,„c w J, 


Kg?] - rSSM 


where 


*<*> = [5?]' 


and 




(102) 

(103) 

(104) 


On the other hand, if we are using the true arc length for horizontal separation, i.e., s = s, 
we can substitute using (38) to obtain, 

, x . , x [ s/a fHij'( s )-\ 

= sign (<*) sign(^) | [-^\ c w 

s/a fPij'i 8 )] , , * r/ x f/ ,, ( s )]l 1 

[^ML-^7] cos(s/a) - H^H] CiyCyj ) 


sin 2 (s/a) 


= sign (</>,•) sign(^)|fl(s)c yy + T(s)c, y c w }, 


where 


and 


R(s) = 


s/a 

rKf'Wi 

sin (s/a) 

L ~/3 2 s J 


/.ZZ (§\ 

T{S) = sin 2 (V«j^ (S)COS(5/a) ~ [^“D' 


(105) 

(106) 

(107) 


The computational algorithms for all of the various correlations will require the exact 
computation of Sij = (s,j or s tJ ) and then table look-ups for (ji zz (sij), R(sij ), T’(s.j)] 
or [/i”(s,j), R(sij), T(sij)] as appropriate. Note that in Appendix B we showed that 
lim a ^oH zz (s) = limj_^oM zz (^) = 1> and also that l«m,-K>#( s ) = limi-m R(s) = 1- We can 
also show that lim s ^ 0 T(s) = lim^o ?(s) = 2(ac 4 ) 2 r 1 + (ac 2 ) 2 7 + (a(3) 2 r)(l + *) where /3, 7 
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and 7? are defined in (65) - (67). This limit at s = 0 can be used when tabulating T(s) or 
T(s). 




which according to (75) and (76) can be written 


— I5J ^2 + ^»?(1 + c 5 )£> - S7COSC2S + 27C5BCD}. 


Using (73) and (74) this can be written 

2 

= gm 2a 2 c 2 d{d W (l + c 5 ) + 27C5BC} - jim 

and, since Iimi_ f0 [C''/s] = -|c^, we get 

lim o f(J) = 2a 2 c 2 [7 ? (l + c 5 ) + 27 c 5 ]+| a 2 c 2 


= 2(ac 4 ) 2 r/ + (ac 2 ) 2 7 + (a(3) 2 i](l + J). 

o 


(108) 


For the damped cosine parameters given in (58) - (62), we get lim i - +0 T(s) as 219.2. Note, 
however, that in (102) and (105) lim^, c iy c yj = 0 and lim^,- c yy = 1, so that lim = 


In a manner completely analogous to our treatment of we obtain the following formulas 
from (81) - (89) for the background error correlations: 


ZZ 

Vi ii = 


1J 

l™ 

U 


Vij = 


uz 

t*ii - 


zv 

Hi = 


vz 

Hi = 


uu 

Pii 


u™ - 


vu 

Hi = 

vv 

Hj = 


A*”(*ij) ( 109 ) 

-sign (<pj)al 3 R(sij)ciy (110) 

- sign (^,)a/?£(s, ,)<:„• (HI) 

+ sign(^)a/3B(s tj )[^-] (112) 

+ sign(^)a/?B(s,y)[^-] ( 113 ) 

+ sign(v?,^){fl(s,j)c w + f(s t] )c ty c y] } (114) 

-8ig n(Wi ){«( 5j> )[-^] (115) 

-s ign(w) {fl(l,)[^] + f (io )[-±l_] c ,„} ( 116 ) 

+ S ig„( W ){fe J )[_^_] + ffe)[^-] (117) 


C.4 Decoupling of the geostrophic assumption in the tropics 

So far we have made uniform use of the geostrophic assumption as expressed in (80). This 
assumption is clearly not valid near the equator, since it requires a division by / = 2 fi sin ip. 
This leads to the ratio of wind to height errors, a u /a z = g(3/\f\, in (97) becoming unbounded 
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near the equator, and to discontinuities at the equator in the definitions of the correlation 
functions given in (110) - (117). 


In equations (110) - (117) for the height-wind error correlations, the discontinuity of the 
factor sign (9) at the equator is removed, as in Baker et al. (1987), by making the approxi- 


mation 


(a/?) sign (v?) « H(<p) = (a/?) sign(y>){l - exp(-|^|)}. 


(118) 


For the wind-wind correlations in (114) - (117) the factor sign(<pi<^) is replaced with 1, i.e., 
the value it has when the two points are in the same hemisphere. Our final equations for 
the height-height, height-wind, and wind-wind error correlations for the upper-air analysis 
in the tropics (15S < <fi < 15N) are 1 


ZZ 

M.i 

Mf? 

uz 

My 

n ?V 

Mu 


vz 

Mu 


rtj 


Mu 


vv 

MO 


M~(*i) 

~H{<Pj) R(Sij ) c iy 
-H{vi)R{si])c y3 

+H (y>i)A(5o) [ 




COS (fij 


:] 






'yy 


+ ^(^u) c ty c yj} 


^yx 




Mu = 


l COS (pj *• ' ' ' ~ LCOS Pj 


Lcos <Pi COS <Pj 


H— ]} 

J LcOSU?; J * 


COS<Pi^ L COS <Pj 


(119) 

( 120 ) 
( 121 ) 
( 122 ) 

(123) 

(124) 

(125) 

(126) 
(127) 


Note that if we set A, = A* and pi = p* = —90° or 90°, and let (<£>), Aj) — > (p ,Aj) we 
obtain the following limiting values My — 1, My* = My* = My = My — My* = My = 
cos(A* - Ap, My = sign(^) sin(A* - Ap and My = - «gn(y>*) sin(A? - Ap. Contour plots 
of the correlation functions defined in (120) - (127) for A* = 0° and pi € {10 , 40 , 90 } 
are shown a s functions of pj and A j in figures 7 12. 


The wind component error standard deviations given by equation (97) are modified in the 
tropics as follows. In the tropical band, 15°S < p < 15°N, <7^ and <Xj are assigned the 
value given in Table 5. At points in the subtropical bands, 25°S < p < 15 S and 

15° TV < p < 25°N, a t u and aj are defined by 


f u \2 _ (a v ,2 _ [ 25- M l/ )2 ■ [M- 1 5ir g 2 /? 2 K) 2 _\ 
K) - K) "I 10 -\Vtro?) +( 1Q \\ | /25 |2 J’ 


(128) 


where /2s = 2fisin(25°). 

1 Inadvertently, the frozen system that produced the five-year reanalysis did not remove the discontinuity 
across the equator associated with the factor sign(t^) in the wind-wind correlations. This problem has been 
corrected in subsequent versions of the system (see section 9). 
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Figure 8: Graphs of the functions nfj, and /i£ for <pi = 10° and A, 





Figure 9: Graphs of the functions /i**, /x* 5 , /xfj, a nd for </>,• = 40° and A; = 0°. 
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*m> m jf m 


Figure 11: Graphs of the 



is /if/, /if/, /if/, and /if/ for y>, = 90° and A, = 0°. 
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Figure 12: Graphs of the functions /!«*, /i£\ and /*w for y>,- 


90° and 


A, 


V • 
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Appendix D: List of Acronyms 


ACARS 

AG CM 

AIDS 

ARINC 

ASDAR 

CFL 

DAAC 

DAO 

DAS 

EOS 

GCM 

GEOS-1 DAS 

GEOS-1 GCM 

GLA 

HIRS 

IAU 

LAPACK 

MLS 

MSU 

NAVA IDS 
NESDIS 

NMC 

NOAA 

NRC 

01 

PBL 

RAS 

RMS 

SSM/I 

SSU 

STRATAN 
TIROS 
TOVS 
TOVS A 
TOVS B 
TOVS C 
UARS 
UTC 


ARINC Communications Addressing and Reporting System 

Atmospheric General Circulation Model 

Aircraft Integrated Data System 

Aeronautical Radio, Incorporated 

Aircraft to Satellite Data Relay 

Courant-Friedrichs-Levy (commonly used computational 
stability condition) 

Distributed Active Archive Center 
Data Assimilation Office 
Data Assimilation System 
Earth Observing System 
General Circulation Model 

Goddard EOS- Version 1 Data Assimilation System 

Goddard EOS- Version 1 General Circulation Model 

Goddard Laboratory for Atmospheres 

High Resolution Infrared Sounder 

Incremental Analysis Update 

Linear Algebra Package 

Microwave Limb Sounder 

Microwave Sounding Unit 

Navigation Aids (in Table 1 this refers to a special class of pilot balloons) 

National Environmental Satellite, Data, and Information 

Service 

National Meteorological Center 

National Oceanic and Atmospheric Administration 

National Research Council 

Optimal Interpolation 

Planetary Boundary Layer 

Relaxed Arakawa-Schubert 

Root Mean Square 

Special Sensor Microwave Imager 

Stratospheric Sounding Unit 

Stratospheric Analysis system 

Television Infrared Observing Satellite 

TIROS Operational Vertical Sounder 

TOVS clear sky retrieval 

TOVS partly cloudy sky retrieval 

TOVS cloudy sky retrieval 

Upper Atmospheric Research Satellite 

Universal Time Coordinated 
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